/*	This program produces Appendix Figure D.4 */

***** Set directories 
local dir_clean 	"~/Dropbox/Retirement gaming/clean"
local dir_output 	"~/Dropbox/Retirement gaming/output/dataverse"

use "`dir_clean'/mainsample_medbcw.dta", clear


local listcontrolsDD  = " i.t	i.t#i.self_empl i.ndep_cat_1stobs##i.self_empl i.ciiu2_1stobs##i.self_empl i.ndep_cat##i.self_empl 	i.year#i.ciiu2_1stobs i.year#i.ndep_cat_1stobs i.year#i.ciiu2_1stobs#i.self_empl i.year#i.ndep_cat_1stobs#i.self_empl" 
local listcontrolsDDE = " i.t	i.t#i.small_1stobs i.ndep_cat_1stobs##i.small_1stobs i.ciiu2_1stobs##i.small_1stobs i.ndep_cat##i.small_1stobs i.year#i.ciiu2_1stobs i.year#i.ndep_cat_1stobs i.year#i.ciiu2_1stobs#i.small_1stobs i.year#i.ndep_cat_1stobs#i.small_1stobs " 


*** DID SELF-EMPLOYED VS. EMPLOYED LARGE FIRMS  
preserve
replace time_bcw=-5 if time_bcw==-6
keep if (empl==1 &  small_1stobs==0)| (self_empl==1)
*Normalize wages for each group
sum W if time_bcw==-1  &  empl==1
local meanw=r(mean)
g Wz=W/`meanw' 
sum W if time_bcw==-1  &  self_empl==1
local meanw=r(mean) 
replace Wz=W/`meanw' if self_empl==1
bysort self_empl: sum Wz if time_bcw==-1
* Age dummies interacted with self_employed
tab time_bcw, gen(time_bcw_dums)
forvalues d = 1(1)13 {
	g did_edad`d' = time_bcw_dums`d'* self_empl
}		
drop time_bcw_dums5 // drop year before ref age, so everything becomes relative to that
drop did_edad5	
* Regression
reghdfe Wz  time_bcw_dums* did_edad* , absorb(`listcontrolsDD') vce(cluster i )
* Sample size
global n=e(N)
global N=e(N_clust1)
* Coefficients and sd
for any beta sd: gen X = .
forvalues d = 1(1)13 {
	capture qui replace beta = _b[did_edad`d'] if time_bcw==`d'-6
	capture qui replace sd   = _se[did_edad`d'] if time_bcw==`d'-6
} 
replace beta=0 if time_bcw==-1 
replace sd=0 if time_bcw==-1  
* Plot the results
collapse beta* sd*, by(time_bcw)
gen sd_top = beta + 1.96*sd	
gen sd_bot = beta - 1.96*sd
gen sd_top10 = beta + 1.645*sd	
gen sd_bot10 = beta - 1.645*sd
twoway rcap sd_top sd_bot time_bcw, ///
	cmissing(n) lwidth(thin) lcolor(navy)  xsc(r(-5 7)) xlabel(-5(1)7)  xline(-1, lc(green))   ysc(r(-0.03 .09)) ylabel(-0.03(0.03).09,grid)   ///
	|| rcap sd_top10 sd_bot10 time_bcw, cmissing(n) lwidth(thin) lcolor(navy)  xsc(r(-5 7)) ///
	|| scatter beta time_bcw, mcolor(navy) scheme(s1color) lcolor(navy) lpattern(solid) yline(0, lcolor(gs2)) ///
	title("Reported earnings, relative to time -1")  note("N=${n}, Individuals=${N}") ///
	legend(off) ytitle("Estimated Coefficients") xtitle("Years Relative to Predicted Start of BCW")   
graph export "`dir_output'/figureD4a.png",  replace 	
restore

		
 
*** DID SMALL FIRM-EMPLOYED VS. EMPLOYED LARGE FIRMS 
preserve
replace time_bcw=-5 if time_bcw==-6
keep if (empl==1 &  small_1stobs==0 )| (empl==1 &  small_1stobs==1)
*Normalize wages for each group
sum W if time_bcw==-1  &  empl==1 &  small_1stobs==0
local meanw=r(mean)
g Wz=W/`meanw' 
sum W if time_bcw==-1  &  empl==1 &  small_1stobs==1
local meanw=r(mean) 
replace Wz=W/`meanw' if  empl==1 &  small_1stobs==1
bysort small_1stobs: sum Wz if time_bcw==-1
* Age dummies interacted with self_employed
tab time_bcw, gen(time_bcw_dums)
forvalues d = 1(1)13 {
	g did_edad`d' = time_bcw_dums`d'* small_1stobs
}		
drop time_bcw_dums5 // drop year before ref age, so everything becomes relative to that
drop did_edad5	
* Regression
reghdfe Wz  time_bcw_dums* did_edad* , absorb(`listcontrolsDDE') vce(cluster i )
* Sample size
global n=e(N)
global N=e(N_clust1)
* Coefficients and sd
for any beta sd: gen X = .
forvalues d = 1(1)13 {
	capture qui replace beta = _b[did_edad`d'] if time_bcw==`d'-6
	capture qui replace sd   = _se[did_edad`d'] if time_bcw==`d'-6
} 
replace beta=0 if time_bcw==-1 
replace sd=0 if time_bcw==-1  
* Plot the results
collapse beta* sd*, by(time_bcw)
gen sd_top = beta + 1.96*sd	
gen sd_bot = beta - 1.96*sd
gen sd_top10 = beta + 1.645*sd	
gen sd_bot10 = beta - 1.645*sd
twoway rcap sd_top sd_bot time_bcw, ///
	cmissing(n) lwidth(thin) lcolor(navy) xsc(r(-5 7)) xlabel(-5(1)7)  xline(-1, lc(green))  ysc(r(-0.03 .09)) ylabel(-0.03(0.03).09,grid)   ///
	|| rcap sd_top10 sd_bot10 time_bcw, cmissing(n) lwidth(thin) lcolor(navy) xsc(r(-5 7)) ///
	|| scatter beta time_bcw, mcolor(navy) scheme(s1color) lcolor(navy) lpattern(solid) yline(0, lcolor(gs2)) ///
	title("Reported earnings, relative to time -1")  note("N=${n}, Individuals=${N}") ///
	legend(off) ytitle("Estimated Coefficients") xtitle("Years Relative to Predicted Start of BCW")   
graph export "`dir_output'/figureD4b.png",  replace 	
restore


clear all
exit
